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INTRODUCTION 


-^Certain  u;*rs  of  the  Global  Positioning  System  (GPS)  use  the  geodetic  coordinate  system  to  define 
their  receiver  position.  Satellit-  positions,  computed  from  either  ephemeris  or  almanac  data,  are  best 
stated,  however,  in  terms  of  Earth-centered,  Earth-fixed  (ECEF)  cartesian  coordinates.  This  presents  a 
problem  when  computing  values  for  a  dilution  of  precision^jn  position-A-the  measure  by  which  the  four 
satellites  most  suitable  for  range  measurements  are  determined/^  /{x>3  , 


Before  receiving  transmissions  from  GPS  satellites  and  using  ther'  to  determine  precise  receiver 
position,  the  user  must  determine  which  four  satellites  are  in  the  configuration  that  will  maximize  system 
accuracy.  This  is  accomplished  by  computing  a^JUution  of  precision) value  for  all  groups  of  four  satellites 
that  are  in  view.  For  GPS  users  who  measure  their  receiver  position  in  geodetic  coordinates,  this  value 
should  reflect  the  use  of  these  coordinates,  rather  than  ECEF  cartesian  coordinates. 

This  report  outlines  five  algorithms,  all  used  in  the  process  of  computing  ablution  of  precisiop-  in 
Position.  The  method  to  compute  satellite  position  is  described  assuming  the  user  is  using  almanac  data, 
and  the  dilution  of  precision  in  position  term  is  defined  in  terms  of  geodetic  coordinates.  f'-  p  A  ) 

The  basic  procedure,  as  outlined  by  the  algoiithms,  is  as  follows 

Use  almanac  data  to  find  satellite  position  at  time  (t) 

Convert  receiver  position  estimate  from  geodetic  coordinates  to  ECE  -  cartesian  coordinates 

Establish  whether  a  satellite  is  healthy  and  in  view 

Compute  covariance  matrixes  for  all  groups  of  four  satellites  that  are  healthy  and  in  view 

Rotate  each  covariance  matrix  to  be  in  terms  of  geodetic  coordinates,  and  compute  a  dilution  of 

ptocision  in  position 


The  group  of  four  satellites  that  yields  the  smallest  dilution  of  precision  in  position  is  the  group  to  use 
when  making  range  measurements  from  satellite  transmissions. 
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Each  of  the  five  algorithms  details  how  to  accomplish  one  of  the  tasks  previously  stated.  The  it.  .ns 
that  are  needed  as  input  to  the  routines  are  listed  and  defined  first.  Th?  «teps  and  equations  required  to 
produce  the  desired  output  follows.  The  algorithms  are  presented  in  a  straightforward  manner,  and 
computer  code  could  be  written  from  them.* 

The  appendixes  provide  a  more  thorough  explanation  of  the  equations  \nd  procedures  detailed  in  the 
five  algorithms.  The  collection  of  material  found  in  the  appendixes  should  prove  to  be  as  useful  as  the 
algorithms  themselves.*  * 

None  of  the  material  in  this  report  is  developed  for  the  first  time  here,  but  as  far  as  it  is  known,  it  has 
never  been  assembled  as  in  this  report.  The  usefulness  of  having  this  inform-*  i'.oa  combined  into  one 
document  was  a  major  reason  for  assembling  this  report. 


PROCEDURES 

FINDING  SATELLITE  POSITION  FROM  ALMANAC  DATA1 
The  user  obtains  the  following  input  data  from  the  almanac 
ISAT  Satellite  tracker  number 

WNA  Reference  week  of  almanac  data 

e  Eccentricity 

tQa  Reference  time  of  almanac  data  (s) 

iQ  Inclination  angle  (0.30  semicircles) 

6j  Correction  to  inclination  (sc) 

H  Rate  of  right  ascension  (sc/s) 

(ae)!//2  Square  root  of  semimajor  axis  length  (m) 

*  These  algorithms,  or  sections  of  these  algorithms,  alreadx  exist  as  code  or  in  program-design  documentation  These  were 
assembled  from  Reference  I. 

•  *  Nearly  all  the  derivations  collected  in  the  appendixes  can  be  found  in  either  Reference  I  or  2. 
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SlQ  Right  ascension  (sc) 

vj  Arbument  of  perigee  (sc) 

M0  Mean  anomoly  (sc) 

The  user  provides  the  following  input  data 

t  Time  at  which  position  is  to  be  known  (s) 

WN  Week  number  at  which  position  is  to  be  known 
ptg  WGS-72  value  for  the  Earth’s  gravitational  constant 

S2e  WC»S  72  value  for  the  Earth’s  rotation  rate 

The  user  receives  the  following  output  data 
xs*  ys’  zs  Satellite  position 

E  Eccentric  anomoly 

Define  semimajor  axis  length,  ae.  mean  motion,  N,  time  difference  (positive  or  negative)  between  t  and 
toa.  At,  and  a  new  mean  anomoly.  M.* 

ae  =  <ae'/2)2 

N,  =No=0ig/ae3^ 

At  =  t  +  604,800  (\»  N-WNA)  -  tQa 
M  =  E  =  Mq  +  N,  •  At 

Solve  for  eccentric  anomoly  by  repeating  the  next  two  lines  eight  times  or  until  tau  (r)  becomes  lessthanor 
equal  to  1  x  10”10. 

t  =  E  -  e  sin  E  -  M 

E  =  M  +  e  sin  E 

*  Sec  Appendix  A  for  derivations  of  some  of  these  equations.  This  routine  is  part  of  Reference  1:  it  is  3lso  contained  in 
Reference  3. 
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Solve  for  true  anomoly,  v,  corrected  radius,  r,  and  corrected  inclination  i. 
sin  v  =  (1  -  e2  f2  sin  E  /  (1  -  e  cos  E ) 
cosv  =  (cos  E  -  e)/ (1- e  cos  E) 
r  =  ae  (1  -  e  cos  E) 
i  =i0  +  6i 

Find  the  corrected  longitude  of  the  ascending  node  and  the  position  of  the  satellite  in  the  orbital  plane. 
(Let  u  =  f  +  gj) 

SI  =  f20  +  SI  •  At  -  S2e  (At  +  toa) 

xs  =  r  -cosu  =  r  -  (cost'  cos  to  -  sin  v  sin  col 

y '  =  r  •  sin  u  =  r  ♦  [sin  v  cos  u  +  cos  v  sin  gj] 

Solve  for  position  in  the  ECEF  cartesian  frame. 

xs  =  x'$  cos  SI  -  y  s  cos  i  sin  SI 
ys  =  x  s  sin  SI  +  >  s  cos  i  cos  SI 

h  =y/ssini 

CONVERTING  RECEIVER  POSITION  FROM  GEODETIC  COORDINATES  TO  ECEF  CARTESIAN 
COORDINATES1 

The  user  inputs  the  following  data 

X,  h  Geodetic  latitude,  longitude,  and  height  of  receiver 

ae  Earth’s  semimajor  axis  length 

e  Earth's  eccentricity 

The  user  obtains  output  data  in  the  form  of  the  x.  y,  and  z  receiver  coordinates  (figure  1). 
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FIGURE  1.  RECEIVER  POSITION 


If  latitude  and  longitude  are  input  in  degrees,  convert  to  radians. 
<t>  =  0  •  jt  /  180 

X  =  X  •»/  180 

Convert  to  ECEF  cartesian  coordinates. 


x 


y 


z 


■[ 

■[ 

■[ 


Vi -  e2  sin2  <t> 

h, 

Vi  -  e2  sin2  <f> 
(i  -  e2 )  3p 

Vi  -  e2  sin2  $ 


cos  0 cosX 


cos0  sin  X 


+  h 


rin  4> 


DETERMINING  SATELLITE  QUADRANT  AND  ESTABLISHING  WHETHER  A  SATELLITE  IS 
HEALTHY  AND  IN  VIEW1 


The  user  inputs  the  following  data 
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x$,  ys,  zs  Satellite  position  in  ECEF  cartesian  coordinates 

x,  y,  z  Receiver  position  in  ECEF  cartesian  coordinates 

0cut  One  cutoff  angle  for  each  of  the  four  quadrants 

The  receiver  is  used  as  a  reference  point  to  define  the  quadrants  as  follows  (see  Appendix  C  for 
further  explanation) 

0  <  quadrant  1  <  90  North-East 

90  <  quadrant  2  <  180  North-West 

1 80  <  quadrant  3  <  270  South-West 

270  <  quadrant  4  <  360  South-East 

Establish  the  satellite  quadrant,  and  define  the  vector  from  the  receiver  to  the  satellite.  P.  by  its  three 
components  Px.  Pv  and  Pz. 

Px  =  xs  “  x 

Py  =  ys  -  y 
pz  =zs-7 

Define  vectors  E  and  N,  which  with  zenith.  Z,  define  a  new  coordinate  system  centered  at  the  receiver. 
(Ex,  Ey,  E,)-  (-y,  x,  o)  East 
(Nx.  Ny,  Nz)  =  (-xy,  -yz,  xy  +  y2 )  North 

Define  the  magnitudes. 

1  E  I  =  (Ex2  +  Ey2  +  Ez2  V2 
I  N  I  =  (Nx2  +  Ny2  +  Nz2  ft 
Compute  satellite  position  in  the  new  E  — N  plane. 
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PE  =  (-Px  -y.Py  -x,o) 

P  ‘  N  =  (-Px  •  x  •  z,  -Py  •  y  •  z,  Pz  •  x  •  y  +  Pz  •  y2 ) 


(Pf,  PN) 


P^  and  P^  determine  the  satellite  quadrant. 

If  Pg  =0  and  P^  =  0  (exactly  on  azimuth,  Z)  satellite  is  in  quadrant  1 
If  Pg  >0  and  P^  >  0  satellite  is  in  quadrant  1 

If  Pg  <0  and  P^  >  0  satellite  «s  in  quadrant  2 

If  Pg  <0  and  P^  <  0  satellite  is  in  quadrant  3 

If  Pg  >0  and  P^j  <  0  satellite  is  in  quadrant  4 


Establish  whether  or  not  the  satellite  is  in  view. 

(Px.  Py.  P7.)  =  (x$  -  x,  ys  -  y,  zs  -  z) 

Compute  the  following,  where  R  is  the  receiver  position  vector. 

IP  I  =  (Px2  +  Py2  +  Pz2  p 

|  R  |  =  (x2  +  y2  +  z2 

P  •  R  =  PX  •  X  +  Py  *  i  +  P7  ♦  Z 

The  angie  between  P  and  the  horizon  is  given  by  p 


If  sin  p  <  sin  flcut  the  satellite  is  out  of  view 
If  sin  P  >  sin  0cut  the  satellite  is  in  view 
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COMPUTING  THE  COVARIANCE  MATRIXES  FOR  ALL  GROUPS  OF  FOUR  SATELLITES  THAT 
ARE  HEALTHY  AND  IN  VIEW  1 

The  user  inputs  the  following  data 

xs*  ys*  zs  Satellite  positions  for  four  satellites  (and  tracker  number  subscript,  s) 
x,  y,  z  Receiver  position 

For  each  of  the  four  satellites  compute  a  4  x  4  matrix  Ms. 

First  let 

Ax$  =  (xs-  x) 

Ays  =  (ys-  y) 

Azs  =  <zs— z) 

Rs  =  [Ax$:  +  Ays2  +  Azs:  ]'/z 
and 

Pls  =  Axs/  Rs 

p~s  =  *y$!  Rs 

P3S  =  Azs/  Rs 
P4S  =  1 

Now, 


P,s  Pis 

pisp::s 

pls-p3s 

P1S  ' P4: 

p2s-p's 

p2s-p-s 

P2s-p3s 

P“S  ■  P4 

p3s-pl$ 

p?sp-s 

p3s -p3s 

P3S  •  P4 

_p4s-pis 

P4S  •  P2$ 

P4S-P3S 

P4S  •  P4 
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Let  Mtovaj  be  the  sum  of  the  four  M  s 
4 


s=  1 


The  final  step  in  this  procedure  is  to  compute  the  inverse  of  Mtotaj.  This  is  the  covariance  matrix.  In  the 
next  section  this  is  denoted  by  M-1 . 


COMPUTING  A  DILUTION  OF  PRECISION  IN  POSITION  FROM  THE  COVARIANCE  MATRIX 
IN  TERMS  OF  LATITUDE  AND  LONGITUDE  1 


The  user  inputs  the  following  data 


M 


-i 


Covariance  matrix 


<f>,  X,  h  Geodetic  receiver  position 


Eccentricity 


Semimajor  axis  length 


x,  y,  z  Cartesian  receiver  position 


The  user  receives  the  following  data 


DOP^h  Dilution  of  precision  in  position, 


DOP 


$X 


Dilution  of  precision  in  position. 


Let  M_l  be  written  as  follows 


m  1 1  m  1 2  m  i  3  m  1  4 


m 


2  l  m2  2  m2  3  m24 


reflecting  latitude,  longitude,  and  height 
reflecting  latitude  and  longitude 


M-'  = 


m3.m32,n33m34 


m4  I  m42  m43  m44 
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M~*  was  generated  from  variables  in  the  ECEF  cartesian  (x.  y.  z)  coordinate  system.  For  users  measuring 
position  in  this  coordinate  system,  geometric  dilution  of  precision  is  computed  by 

GDOP  =  Tr[M_1]  =  m(  %  +  mJ2  +  mJ3  +  m44 

The  procedure  described  following  these  notes  yields  a  dilution  of  precision  suitable  for  users  measuring 
position  in  the  geodetic  ( 0 ,  X,  h)  coordinate  system.  Two  dilution  of  precision  equations  are  defined;  one 
for  users  always  at  nearly  constant  height  (DOP^  ),  and  one  for  users  whose  height  varies  (DOP  ^j,)- 

Form  a  new  matrix  by  "rotating”  the  upper  left  3x3  matrix  in  M_1  into  the  geodetic  coordinate  system. 
Equivalently,  multiply  the  upper  left  3x3  matrix  in  M'1  by  the  following  matrix  of  partials 


30 

30 

30 

3x 

3y 

3z 

3X 

3X 

3X 

3x 

3y 

3z 

3h 

3h 

3h 

3x 

ay 

3z 

The  following  procedure  shows  the  equations  for  the  partials  in  the  above  matrix.  Appendix  D  details  how 
they  are  derived. 


Compute  the  range  in  the  x— y  plane  and  the  radius  of  curvature. 
R  =  Vx-  +  y2 


RADC 


ae  (i  -  e2 ) 

[ i  —  e2  sin2  0l  ^ 


+  h 


Comp.ue  partials  of  geodetic  latitude. 

30  _  -cos  X  sin  0  30 

3x  RADC  3y 


-sin  X  sin  0 
RADC 


Compute  partials  of  longitude. 

3X  _  -sinX  3X  _  cosX 

3x  R  3y  R 


30  cos  0 


3z  RADC 


3X 

- —  =  zero 
3z 
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Compute  partials  of  height, 

ah 


3x 


=  cos  X  cos  0 


Dilution  of  precision. 

dfn. 


DOP 


0\ 


DOP 


0Ah 


3h 

—  =  sin  X  cos  0 

ay 


—  + 


3x 


DOP^X  + 


ah 

—  =  sin  0 

3z  v 


dV. 

30mu 

dX: 

3H: 

1 

ay 

37 

3x 

ay 

-r 

az 

av, 

aS, 

i 

3h„ 

1 

ax 

dy 

37 
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APPENDIX  A 

USING  ALMANAC  DATA  TO  COMPUTE  SATELLITE  POSITION 


A-l 
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ORBITAL 

PATH 


i  =  INCLINATION 
w  =  ARGUMENT  OF  PERIGEE 
n  =  LONGITUDE  OF  ASCENDING 
NODE 


FIGURE  A-l.  CLASSICAL  ORBITAL  ELEMENTS 


When  computing  satellite  position  from  almanac  data,  first  find  the  values  of  the  sine  and  cosine  of  the 
true  anomoly.'-1 

cos  v  =  (cos  E  —  e)/(  1  -  e  cos  E) 
and 

sin  v  =  ^/l  —  e 2  •  sin  E/(  1  —  e  cos  E) 

True  anomoly  is  the  angle  shown  in  Figure  A-2. 


''■iR.  R.  Bate.  D.  D.  Mueller,  and  J.  l-‘.  While.  f'n>nl<nnrnnil\  »/  As.rndwiamii'-  (No'  York'  Dover  Publications.  Inc..  I9'I). 
pp.  183-187. 
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FIGURE  A-2.  TRUEANOMOLY 


Note  two  things 

1 .  For  any  ellipse,  like  the  one  in  Figure  A-2 .  the  focus  is  positioned  such  that 
b2  =  a2  -  c2 

b2  =  z2  -  e2a2 

and 

b  =  a  (i-  c2  f'2 

2.  Given  two  points  with  the  same  x  coordinate,  one  on  the  circle  and  one  on  the  ellipse  as  drawn 
in  Figure  A-3,  the  ratio  of  y  compcnents  is 

y ellipse  _  ^ 
y  circle  a 


A-4 
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FIGURE  A-3.  ELLIPSE  INSCRIBED  WITHIN  A  CIRCLE 


This  can  be  seen  from  the  equations  for  an  ellipse  and  a  circle.  For  a  circle 


Using  Figure  A-4,  solve  for  r,  then  solve  for  true  anomoly. 


A-5 
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K 


FIGURE  A-4.  TRUE  ANOMOLY 


r2  =  x2  +  Py2 

r2  =  [c  -  a  cos  El 2  +  (  ^  ^  •  a  sin  El 2 

r2  =  a2  e2  -  2a2  c  cos  E  +  a2  cos2  E  +  b2  sin2  E 

r2  =  a2  [e2  -  2e  cos  E,+  cos2  E  +  ( 1  -  e2  )  sin2  El 

r2  =  a2  [  e2  -  2e  cos  E  -  c:  sm2  E  +  11 

Use  the  expression  ( 1  —  e  cos  E)2  =  1  —  2c  cos  E  +  e2  cos2  E 
r2  =  a2  [e2  +  (1  -  ecosE)2  -  e2  cos2  E  -  c2  sin2  El 

so 

r2  =  a2  ( 1  -  e  cos  E)2 

and 

r  =  a  ( 1  -  e  cos  E) 
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Finally 


cos  v  - 


-  sin  0  =  -  — 
r 


cos  v  -  - 


(c  -  a  cos  E)  _  cos  E  -  e 
a(l  -  e  cos  E)  1  -  e  cos  E 


sin  v  - 


Py  b  sin  E 
cos  0  =  —  =  - 


_  V 1  -  e2 


sin  v  -  — 


1  -  e  cos  E 


Using  Figure  A-5.  define  satellite  position  in  the  orbital  plane. 


"o  =  TRUE  ANOMOLY 

gj  =  ARGUMENT  OF  PERIGEE 


con  tc 


,  ORBITAL 
I  PATH  &  PLANE 


*•  PERIGEE 


FIGURE  A-5.  SATELLITE  POSITION  IN  THE  ORBITAL  PLANE 


Define  a  coordinate  system  in  the  orbital  plane  where  the  x'  axis  lies  alorg  the  ascending  node,  and  the 
y'  axis  is  perpendicular  to  the  x/  axis  with  the  origin  at  the  focus.  See  Figure  A-5. 
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Let 

uk  = 
so  that 

rx<  =  x‘  cos  uk  =  x'  cos  v0  cos  to  -  r'  sin  vQ  sin  to 
Ty/  =  x‘  sin  uk  =  r'  sin  vQ  cos  to  +  x'  sin  to  cos  uQ 

Recall  that 

x'  -  a  (1  -  e  cos  E) 

The  corrected  longitude  of  tne  ascending  node  is  illustrated  by  Figure  A-6. 


FIGURE  A-6.  CORRECTED  LONGITUDE  OF  THE  ASCENDING  NODE 


A-8 
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=  Initial  right  ascension,  measured  from  Greenwich,  at  the  beginning  of  the  week,  to  the 
ascending  node  at  toa 

ft  =  Rate  of  right  ascension,  the  rate  of  orbital-plane  drift  around  a  fixed  Earth 
fte  =  Earth’s  rotation  rate 

t  =  Time  at  which  satellite  position  is  to  be  known 
tQa  =  Time  of  almanac  data 
ftk  ~  Longitude  of  ascending  node 
From  Figure  A-6,  ftk  is  given  by 

ftk  —  ft0  +  ft  (t  —  toa)  —  fte  (t  —  tQa) 
ftk  =  no+  &  (*  ~  loa )-  fte  ' 1 

and  where 

lk  =  t  ~  +  604,800  •  (#  of  weeks  between  t  and  toa) 

^k  “  ^o+  ^  ’  Me  ~  ^e  ^k  +  *oa^ 

Figure  A-7  illustrates  the  conversion  from  position  in  the  orbital  plane  to  position  in  ECEF  cartesian 
coordinates. 

rx  =  rx/cosftk-  ry,cosik  sin  ftk 
ry  =  rx/sin  I2k+ Ty/cosik  cosS2k 

r7  =  Ty/sinik 

where  r  /  and  r  ,  are  the  orbital  plane  coordinates,  as  defined  previously, 
x  y 
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APPENDIX  B 

CONVERTING  GEODETIC  COORDINATES  TO  ECEF  CARTESIAN  COORDINATES 
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Before  converting  geodetic  coordinates  to  ECEF  cartesian  coordinates.8  1  note  two  things 
1.  As  seen  in  Figure  B-l.  for  any  ellipse,  the  focus  is  positioned  such  that 
b2  =  a2  -  c2 
b2  =  a2  —  e2  a2 
where 

e  =  c  /  a 
and 

b  =  a  (1  -  e2^ 


FIGURE  B-l.  ELLIPSE 


B-1r.  R.  Bate.  0.  D  Mudlcr.  and  .1  I-  White,  Fundamental*-  Aerodynamics  (New  > ork  Dover  Publicauons.  Inc...  I9'|). 
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2.  For  a  point  on  a  circle  and  a  point  on  an  inscribed  ellipse  with  the  same  x  component,  as  seen  in 
Figure  B-2,  the  ratio  of  y  components  is 

yellipse  _  ^ 
y circle  a 


FIGURE  B-2.  ELLIPSE  INSCRIBED  WITHIN  A  CIRCLE 


From  Figure  B-3  and  the  notes  on  the  previous  page 
Pj  =  ae  cos  0 

and 


ae  sin  P  =  ac  v/'j  c:  •  sin  ^3 


where  these  are  to  be  redefined  in  terms  of  6.  not  3  . 
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j  =  AXIS  FROM  EARTH’S  CENTER 
THROUGH  THE  NORTH  POLE 
i  =  AXIS  IN  EQUATORIAL  PLANE 
THROUGH  LONGITUDE  OF  P 
*e  =  SEMIMAJOR  AXIS 
be  =  SEMIMINOR  AXIS 
0  =  GEODETIC  LATITUDE 
P  =  GEOCENTRIC  LATITUDE 
7  =  TANGENT  TO  EARTH’S 
SURFACE  AT  P 
n  =  NORMAL  TO  t  FORMS 
GEODETIC  LATITUDE 


FIGURE  B-3.  CROSS  SECTION  OF  EARTH  INSCRIBED  WITHIN  A  CIRCLE 


dj  A  di 

The  slope  of  t  is  —  and  the  slope  of  n  is  -  —  =  tan  d 
di  dj 


so  differentiating  Pj  and  Pj  yields 


Now  find  cos  /?  and  sin  /J  in  terms  of  d  . 
Let 


B-5 


V^WiVTST^^ii-.A"*  ViTW^VtXW.'gt K%1  ’ -~  VL'Cvr’-».i»s_^»tx 
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so  that 


sin  P  = 


sin  0  = 


ard 


y 


A2  +  B2 


y/l  ~  e2  • 


sin  ^ 


yr 


e2  sin2  0 


ae  cos  0 


/l-e2 


sin2  0 


and 


cos  /3  = 


B 


cos  0  = 


v/a2  +  B2 


cos  0 


yr 


e2  sin2  0 


ag  (1  -  e2  )sin  0 


/ 


1  -  e2  sin2  0 


Let  P'  be  a  height  h  above  the  Earth's  surface,  as  in  Figure  3-4.  From  this,  it  is  seen  that  the  coordinates 
of  P '  in  the  i-j  frame  are 


yr 


+  h 


L  v  1  -  e2  sin2  0 


COS0 


aP  ( 1  e2  ) 


I  yrr.. 

L  V  1  -  e*  sin* 


+  h 


sin  0 


Now  place  the  i-j  frame  onto  an  x-v-z  frame,  with  the  j  and  z  axes  coinciding  as  in  Figure  B-5.  The 
ECFF  cartesian  coordinates  are 


Px  -  Pj  cos  X  = 


ae 


L  V  I  e:  sin2  0 


+  h 


cos  0  cos  X 


Py  =  Pj  sin  X  = 


ae 


_  y  1  e2  sin2 


+  h 


cos  0  sin  X 


( 1  -  e2  )a„ 


L  \/  1  e2 


sin*  0 


+  h 


sin  0 


B-b 
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FIGURE  B-5.  ECEF  CARTESIAN  COORDINATES 


B-8 
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APPENDIX  C 

DETERMINING  SATELLITE  QUADRANT  AND  WHETHER  THE  SATELLITE  IS  IN  VIEW 
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Figure  C-l  shows  the  Earth’s  northern  hemisphere  and  the  ECEF  cartesian  coordinate  system.  Receiver 
position  is  given  by  R;  satellite  positior,  is  given  by  S.  (These  data  are  necessary  to  determine  the  satellite 
quadrant  and  whether  the  satellite  is  in  view.r'1 

—a.  ^ 

P  =  S  -  R 


z 


FIGURE  C-l.  NORTHERN  HEMISPHERE  AND  ECEF  CARTESIAN  COORDINATES 


A 

Let  k  be  the  unit  vector  along  the  z  axis.  Define  a  new  cartesian  system,  with  its  origin  at  the  receiver,  by 
taking  the  following  cross  products. 

E  =  k  X  R 
N  =  R  X  E 
Z*EXN 


^  'IS.  i  Mewrhofl.  C ISA R  Formulation  ami  Soiiwarf  Working  Pa/n  r\,  NSWC.  Dahlgreri.  Virginia,  Jan  i 
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E  is  tangent  to  the  Earth’s  surface  at  the  receiver  and  points  eastward. 

N  is  nearly  tangent  to  the  surface  at  the  receiver  and  points  northward.  If  the  Earth  were  a  true  sphere,  N 
would  be  exactly  tangent  at  the  receiver. 

Z  is  perpendicular  to  E  and  N,  and  points  toward  the  zenith. 

The  vectors  E,  N.  and  Z  define  the  four  quadrants,  as  shown  in  Figure  C-2.  The  view  is  down  along  the  Z 
axis.  Let  e  and  n  be  unit  vectors  in  the  east  and  north  directions.  Thus,  satellite  position  relative  to  the 
receiver.  P.  is  given  in  this  system  by 

=  I  P  I  cos  0  c  +  |  P  I  cos  a  n 

=  <Plt P  -El  e  +  IP | [ P -N  I  n 
IPIIEI  IPIINI 

=  Life  +  P  •  Nn 

IE  I  INI 


N 


FIGURE  C-2.  ILLUSTRATION  OF  QUADRANTS 
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Define 

penz  ~  <Pe‘PN> 

By  examining  Pg  and  P^,,  (with  predetermined  boundary  conditions),  the  quadrant  in  which  the  satellite 
is  located  can  be  determined.  Each  quadrant  has  a  certain  predetermined  cutoff  angle,  #cut,  measured 
from  the  horizon.  If  the  angle  between  the  horizon  and  the  vector  from  the  receiver  to  the  satellite  is  less 
than  0  t,  the  receiver  is  deemed  to  be  out  of  view. 

To  determine  whether  the  satellite  is  in  view,  find  0,  as  illustrated  in  Figure  C-3,  and  see  if  it  is  less  than 
or  greater  than  0cut- 


FIGURE  C-3.  TWO  VIEWS  OF  SATELLITE  POSITION 


C-5 
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Where  R  and  P  are  both  known 


sin  j3  =  cos  a 


R  P 
I  R  I  IP  I 


If 


sin  (5  <  sin  0CI.(  the  satellite  is  out  of  view. 
If 


sin  /?  >  sin  0cut  the  satellite  is  in  view. 


C-6 
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APPENDIX  D 

DIFFERENTIATING  EXPRESSIONS  FOR  LATITUDE,  LONGITUDE,  AND  HEIGHT 

WITH  RESPECT  TO  X,  Y,  AND  Z 


NSWC  TR  85-151 


Finding  expressions  for  geodetic  latitude,  longitude,  and  height;  and  differentiating  them  with  respect  to 
x.  y.  and  z  requires  finding  an  expression  for  tan  <j>. 15-1 

As  shown  in  Appendix  B 

f  (  1  -  e:  )  ae  ~| 


z 


I  -  e2  sin2 


sin  <t> , 
1 


and  as  shown  in  Appendix  F 


tan  d, 
tan  d,  = - “ 

I  -  e2 

Ustnu  Figures  D-l.  D-2.  and  the  relations  that  follow.,  find  expressions  for  .  and  . 

dx  dy  dz 


Az' 


f 


i' 


* 


ae  =  SEMIMAJOR  AXIS 
h  =  HEIGHT  ABOVE  EARTH  S 
SURFACE 

=  GEODETIC  LATITUDE 
<t>7  =  GEOCENTRIC  LATITUDE 


FIGURE  D-l .  GEODETIC  AND  GEOCENTRIC  LATITUDES 


I*  *S  t  Mf.nh.'tt  iitSAR  iinJ  S W'nrkirv  Puptrs.  NSWC.  Dahlgrcn.  Virginia.  Jan  1983. 
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Az'cos  =h  sin  (dt  *dj) 


FIGURE  D-2.  HEIGHT  AND  GEODETIC  AND  GEOCENTRIC  LATITUDES 


tan  vJ2  z'  1  z  -  A  z  1 

tan  Q  =  - =  —  •  -  =  -  •  - 

1  -  e2  R  (1  -  e2  )  R  (1-e2) 


and 


h  sin  (d,  -  <fi2) 

Az'  = - =  h  sin  d,  -  h  tan  d2  cos  d, 

cosd2 

Az'  =  h  sin  d,  -  h  •  U  -  e2  )sin  d, 

Az'  =  e2  -h  •  sin  <j>l 


tan  d, 


tand,  - 


Az' 


R  (1  -  e2 )  R  ( 1  -  e2 ) 


+  h 


L  \  1  —  e2  sin2 


1  -  e2  sin2  d, 


Ryrr 


sind, 


R 


ae  sin  d,  h  sin  d,  h  e2  sin  d 

—  +  - - 


e2  sin2  d,  R  (1  -  e2 )  R  <  1  —  e2  > 


D-4 
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Note  that  an  equivalent  expression  is 


a*.  e:  sin  0, 

z  _ 1 

tan  0  =  -  +  / - 

R  R^  1  -  sin2  0, 


From  this  expression.  — -  .  — .  and  can  be  found  as  shown  on  the  following  pages. 
3x  3y  3z 


The  above  equation  for  latitude  has  the  problem  of  containing  0,.  on  both  sides.  To  find  a  numerical  value 
for  0,.  when  z.  e.  a  and  R  =  (x:  +  y2)Vl  are  known,  first  estimate  6,.  by 


0,  -  tan-  (|) 

Insert  this  into  the  right  hand  side,  reevaluate  tan  0(.  and  repeat  until  negligible  change  occurs  between 
iterations. 

30 

Find  — - 
3x 

Use  the  following 


v/x:  +  y: 


ae 


\  /  1  ^ 
v  1  e*  sm* 


+  h 


COS  0 


-  (L)‘-x- 
3x  '  R  >  R3 


x  3  »  ,  e-sin0|cos0i 

-  ( !  e:  sm:  0,  )  :  = 


3x 


( 1  e:  sm-  0(  ) 


/: 


Factor  1/R  out  of  the  expression  for  tan  0( .  and  solve  from  there 


d0, 

d\ 


tan  0( 


tan  0( 


7 

R 

I 

R 


a0  e:  sm  0( 


Rs j 


!  e:  sm:  0, 


ac  •  e*  sin  0( 


7  + 


./V  ,  , 

*  I  e*  sm*  o 


i 
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Differentiate 


r  ..  .2 


“1 


r 


*4  2 
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-x  tan  0, 


-  x  tan  o, 


ae 

ae 

/  . ♦  “ 

V  1  -  e2  sin2  0, 

V^l  -  e2  sin2  0, 

t1  COS2  0, 

1  -  - 

1  -  e2  sin2  0, 


+  h 


dO, 

dx 


ac 

ae 

1  -  e2 

/  +  n 

V  1  -  e2  sin2  0, 

n/  1  -  e2  sin2  o, 

1  -  e2  sin2  0, 

do, 

dx 


Since  x  =  R  cos  X  and  R  — 


/  ; 

v  1  -  e2  sin2  0, 


+  h 


cosO, 


and  letting  the  radius  of  curvature  be 
(1  e2)ae 

RADC  =  /  +  h 

v  1  -  c2  sin2  0, 


It  follows  that 

-xtanO,  =  -RcosXtano, 


-  R  cos  X  tan  0, 


dO 


-  •  RADC • 

cosO, 

-cos  X  sin  0, 


do 

dx 


dx 


RADC 
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Similarly,  for  — ,  the  last  steps  are 

ay 


-ytan^j  =  -R  sin  X  tan  ^ 


-R  sin  X  tan  ^  = 


COS0J 


3«, 

RADC  — 

3y 


sin  X  sin  0j 
RADC 


Now  find  — ,  (let  0  =  0,) 
dz 


tan0  =  -  z  + 
R 


as  before,  but  R  = 


ae  •  •  sin  0 


-  e2  sin*  0 


x2  +  y3 


ae 

=  . - - —  *  i. 

_  v  1  -  e2  sin2  0 


+  h  cosp 


is  now  constant 


Differentiate 


sec2  0 


sec2  0  — 


dr  ac -e2  sin2  0cos0 
—  ^  ■  ■  ■■  •  ■  1  + 
dz  (1  -  e2  sin2  p)3/2  dz 


a^  e2  cos  0 


1  -  e2  sin2  0 


ae  e2  cos  0 


e2  sin2  0 


—  +  1 


1  -  e2  sin2  0  I  (1  -  e2  sin2  0) 


*c2  0  — 


,  r  ^ 

-  ii+  / — 


1  -  e2  sin2  0  (1  -  c2  sin2  0) 


30 

dr 


D-8 
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"J \  -  c1  sin2 


+  h 


0 


cos*  30 


ag  e2  cos  0 


=  1  + 


COS2  0 


37 


J 1  -  e 


2  sin2  0 


30  30 

—  +  h  —  =  cos0  + 
3z  dz 


yrr 


V  1  -  e*  sir 


1  -  e2  sin2  0 
<1  e2  )a. 


e2  cos*  0 
1  -  e2  sm2  0 
1 


30  30 

—  +  h  —  =  cos0 
37  37 


+  h 


30 


cos  0 


Thus 


V 

1  -  e2  sin2  6 

( 1  -  e2  sin2  0) 

J  37 

30 

COS0 

3z 

RADC 

where  again  RADC  is  the  radius  of  curvature 
(1  e2  )  ag 


RADC  = 


( 1  e:  ^n2  0) 


+  h 


3X  3X  3X 

Find  — .  — ,  and  — 

3x  3y  dz 

Since  R  and  the  angle  X  swept  out  by  R  are  in  the  x-y  plane. 
3X 


=  0 


37 


3X  3X 

Now  find  —  and  — 

3x  3y 


x 

R 


cos  X  =  -  =  (  x  '  2  ) ' 2  •  (x2  *  y2  )  =  ( 1  ♦  >  2  x' 2  ) 


d  0 


y/ \  -  e2  sin2  0  (1  -  e2  sin2  0) 

ag  e2  cos2  0  1 


1  -  e2  sin2  0  (1  -  e2  sin2  0) 


d0 

dz 


5 

1 

s 

•V 

S: 

>. 

y. 

« 


,N 

A 


,N 

A 


V 

V 

< 
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s 
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After  differentiating: 


-sin  X 

3X 

3x 

1 

=  - (1  +  y2  x-2)  (-2x“3  y2) 

2 

-sin  X 

ax 

3x 

=  (1  +  y2  x'2)3/2  (xj)3/j  y2 

-sin  X 

ax 

ax 

=  (x2  +  y2)3/2y2 

-sin  X 

ax 

dx 

y: 

R3 

ax 

R  y2  y 

ax 

y  R3  R2 

ax 

sin  X 

— 

=  _ - 

ax 

R 
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UseR  = 


+  h 


2  cm2 


sm*  $ 


30 

cos  0  —  = 

dx 


-cos  X  sin  0 


RADC 


x  = 


and  RADC  = 


(1  -  e1)^ 


(1  -  e2  sin2  0)3^2  +  h 


Start  with 


h  -  / - +  - 

V  1  -  e2  sin2  0  cos* 


and  differentiate 


dh 

-ae  * 

e2  sin  0  cos  0  b0  x 

a 

1 

R  sin  0  d0 

3x 

(1  -  e2 

sin2  0) 3^2 

dx  R 

cos0  cos2  0  dx 

dh 

a,.  •  e2  sin2  0  cos  0  cos  X  cos  X 

dx 

RADC 

*  (1  -  e2  sin 

2  0  )/*  COS0 

>/ 1  -  c2  sin2  0 

3h 

cos  X 

ac  e2  sin2 

0  cos2  0 

\  sin2  0 

3x 

COS  0 

RADC  •  (1  -  e2  sin2  0 ) 

h  ' 

'J  1  -  e2  sin2  0 

3h 

cos  X 

sin2  0 

ae 

e2  cos2  0 

dx 

COS  0 

RADC 

7,  .  . 

V  1  c*  sin 

(1  -  e2  sin2  9) 

dh 

cos  X 

sm2  0 

(<?:  i)ac 

1 

3x 

COS  0 

RADC 

I  -  e2  sin2  0 

(1  -  e2  sin2  0) 

dh 

cos  X 

sin2  0 

RADC  +  I 

dx 

COS0 

RADC 

dh 

dx 

= 

COS  X  COS0 

D-ll 


R  cos  X 


sin2  0  cos  X 

h  - 

RADC  -  cos# 


I  h  sin2  0 

- +  j 

RADC  RADC 


1  -  h  +  I 


h  +  1 
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dh 

Similarly  —  =  sin  X  cos  0 

3y 

dh 

Now  find  — .  where  as  before 
dz 


So 


-a„ 


R 


VI  -  e2 


sin2  0 


COS0 


30 

COS0 

dz 

RADC 

dh 

dz 

3h 

dz 

dh^ 

dz 

dh 

dz 


-?e  e2  sin  0  cos  0  30  R  sin  0  30 


(1  -  e2  sin2  0)3'/2  dz  cos2  0  3z 


-af  •  e2  sin  0  cos2  0 
RADC  •(!  -  e2  sin2  0)3/2 


y/ 1  -  e2  sin2 


+  h 


0 


sin  0 
RADC 

sin  0 
RADC 


|_  y/ 1  -  e2  sin2  0 
(1  -  e2 lag 


e2  cos2  0 

j  _ -  |  +  h 

(1  -  e2  sin2  0) 


\/  1  -  e2 


sin2  0 


(1  -  e2  sin2  0) 


J 

sin  0 
RADC 


sin  0 
RADC 


•  RADC 


and  finally 
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APPENDIX  E 

THE  RELATIONSHIP  BETWEEN  GEODETIC  AND  GEOCENTRIC  LATITUDE 


E-l 
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K-l  ...  , 

The  relationship  between  geodetic  and  geocentric  latitude  '  is  shown  in  Figure  E-l.  Assuming  that  the 
Earth  has  the  shape  of  an  ellipsoid,  show  that 

tan  02 

tan  <4  =  - 

(1  -  e2 ) 


j 


FIGURE  E-l.  RELATIONSHIP  BETWEEN  GEODETIC  AND  GEOCENTRIC  LATITUDES 

As  shown  in  Appendix  A,  where  ihe  x  components  are  the  same,  the  y  components  of  a  point  on  a  circle 
and  a  point  on  an  inscribed  ellipse  have  the  following  ratio 

y  ellipse  _  k 
ycircle  a 


H-lR.  R  Bate.  D.  D.  Mueller,  and  J  h  White.  Fundunniuah  <>t  A\irnd\namus  (New  York  Dover  Publieations.  Ine...  1971). 

pp 
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where  b  =  a  ^/l  -  e2 
Thus,  the  angles  <j>2  and  03  are  related  by 
tan  <(>2  =  / 1  -  e2  tan  03 


_  3y 

The  slope  of  tangent  vector  t  is  — , 

3x 


3x 

and  the  slope  of  its  normal  r  is - 

3y 


From  Figure  E-l 


^e 

x  =  ae  cos  y  =  —  •  ag  sin  <j>3 

ae 


Thus 


3x 

tan  <6  -  -  — 

3y 


ae  sin  <p3 

ae  "s/l  -  e2  cos  03 


tan 

\/ 1  -  e2 


and 


tan 

tan  <f>  = - 

1  —  e2 


E-4 
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